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We use the replica method to study the ideal glass transition of a liquid of identical Hard Spheres. 
We obtain estimates of the configurational entropy in the liquid phase, of the Kauzmann packing 
fraction ipx, in the range 0.58 -j- 0.62, and of the random close packing density tp c , in the range 
0.64 -j- 0.67, depending on the approximation we use for the equation of state of the liquid. We also 
compute the pair correlation function in the glassy states (i.e., dense amorphous packings) and we 
find that the mean coordination number at ip c is equal to 6. All these results compare well with 
numerical simulations and with other existing theories. 



I. INTRODUCTION 

The question whether a liquid of identical Hard 
Spheres undergoes a glass transition upon densification is 
still open [1-4]. If crystallization is avoided, one can ac- 
cess the metastable region of the phase diagram above the 
freezing packing fraction <pf = 0.494, where tp = N £y , 
D is the Hard Sphere diameter, N is the number of parti- 
cles and V is the volume of the container. In this region 
the dynamics of the liquid becomes slower and slower 
on increasing the density. The particles are "caged" by 
their neighbors, and the dynamics separates into a fast 
rattling inside the cage and slow rearrangements of the 
cages. The typical time scale of these rearrangements 
increase very fast around ip g ~ 0.56 and many authors 
reported the observation of a glass transition at these 
values of density [5, 6]. 

If the radius of the cages is sufficiently small and if 
the typical time scale of cage rearrangements is suffi- 
ciently large, the system vibrates around configurations 
that are stable for a very large time and can be threated 
as metastable states. It is then natural to separate the 
total entropy of the liquid in a "vibrational" contribution, 
that accounts for the entropy related to the rattling of the 
particles around the metastable structure, and a "con- 
figurational" entropy that is the number of metastable 
states accessible to the liquid at the considered value of 
density [7, 8]. For many simple potentials such as the 
Lennard-Jones [9, 10] and for more realistic systems as 
well [11, 12] the extrapolation of the measured configura- 
tional entropy at higher density (or lower temperature) 
indicates that there exists a density, called Kauzmann 
density <pk, where the configurational entropy vanishes. 
The system freezes in the lowest free-energy states and 
no more rearrangements of the structure are possible. 
This transition is commonly called ideal glass transition 
or Kauzmann transition [8-13]. Note that the Kauzmann 
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density is expected to be larger than the experimental 
glass transition density, as at <pk the relaxation time 
is expected to diverge so that the system freezes in a 
metastable state, on the experimental time scale, for a 
density ip g smaller than tpx. The density tp g where the 
real glass transition happens (weakly) depends on the ex- 
perimentally accessible time scale. Few estimates of the 
configurational entropy for Hard Spheres are currently 
available [3, 14, 15] and indicate a value of (fx in the 
range 0.58 + 0.62. 

A related problem is the study of dense amorphous 
packings of Hard Spheres. Dense amorphous packings 
are relevant in the study of colloidal suspensions, gran- 
ular matter, powders, etc. and have been widely stud- 
ied in the literature [16-23]. The amorphous metastable 
configurations described above provide examples of such 
packings: when the system freezes in one of these states, 
if one is still able to increase the density in order to re- 
duce the size of the cages to zero (for example by shaking 
the container [17, 18] or making use of suitable com- 
puter algorithms [19, 20, 23]), a random close packed 
state is reached. The problem of which is the maximum 
value of density p c that can be reached applying this 
kind of procedures has been tackled using a lot of differ- 
ent techniques, usually finding values of p c in the range 
0.62 -r 0.67. Another interesting problem is to estimate 
the mean coordination number z, i.e. the mean number 
of contacts between a sphere and its neighbors, in the 
random close packed states. Many studies addressed this 
question usually finding values of z ~ 6. 

Recently, the replica method [13, 24, 25] has been suc- 
cessfully applied to the study of the ideal glass transition 
in simple liquids as the Lennard-Jones liquid. Reliable 
estimates of the configurational entropy, of the Kauz- 
mann temperature and of the thermodynamic properties 
of the glass have been obtained from first principles in 
this way [9, 13, 26, 27]. However, for technical reasons 
this approach could not be extended straightforwardly 
to the case of Hard Spheres; indeed at some stage is was 
assumed that the vibrations around the equilibrium po- 
sitions were harmonic in a first approximation. This ap- 
proximation is not bad for soft potentials, but it clearly 



2 



makes no sense for hard spheres. A related but different 
approach was used in [14], obtaining a reasonable esti- 
mate of the Kauzmann density tpx ~ 0.62; however, the 
estimate of the configurational entropy was wrong by two 
orders of magnitude and the thermodynamic properties 
of the glass could not be computed within this approach. 

The aim of this work is to adapt the replica method 
of [13] to the case of the Hard Sphere liquid, and in gen- 
eral of potentials such that the pair distribution function 
g(r) shows discontinuities. This allows us to compute 
from first principles the configurational entropy of the 
liquid as well as the thermodynamic properties of the 
glass and the random close packing density. We find a 
very good estimate of the configurational entropy that 
agrees well with recent numerical simulations [3, 15], a 
Kauzmann density in the range 0.58 -j- 0.62 (depending 
on the equation of state we use to describe the liquid 
state), and a random close packing density in the range 
0.64 -T- 0.67. Moreover, we find that the mean coordina- 
tion number in the amorphous packed states is z = 6 
irrespective of the equation of state we use for the liq- 
uid, in very good agreement with the result of numerical 
simulations [19, 20, 23]. 

The structure of the paper is the following: in section II 
we outline the replica method of [13]; in section III we 
show how it can be adapted to the case of Hard Spheres; 
in section IV we resume the main formulae from which we 
derive our results; in section V we present our main re- 
sults about the configurational entropy of the liquid and 
the thermodynamic properties of the glass; in section VI 
we discuss the behavior of the correlation functions in 
the glass phase; finally, in section VII we compare our 
results with previous works. 



II. THE REPLICA APPROACH TO THE 
STRUCTURAL GLASS TRANSITION 

The replica method was successfully adapted to the 
study of the glass transition of simple liquids in a series 
of recent papers [9, 13, 25-27] . The strategy as well as the 
physics beyond it have been described in detail in [13]: 
in this section we will only review the main steps of this 
approach in order to establish some notations. 



A. The molecular liquid 

Let us consider here a system at fixed density as in [13]. 
The discussion is trivially extended to the case of interest 
here where the density is the control parameter. 

Close to the glass transition the phase space is discon- 
nected in an exponential number of states. The number 
of states of free energy / is called N{f) = cxpVE(/). 
The complexity E(/) is a concave function of / and van- 
ishes at some value f m in- One can write the partition 



function Z in the following way: 

z = e -[3NF(T) _ e -P N f<* 
a 
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where /* is such that /3<I>(/) — f3f — E(/) is minimum. 
The ideal glass transition is met at the temperature Tk 
such that f*(T K ) = f min and £(/*) = 0. 

The basic idea of the replica approach [13, 25] is to 
consider to copies of the original system, constrained to 
be in the same state by a small attractive coupling. The 
partition function of the replicated system is then 



Z m = e -0N^(m,T) ^ ^ e -f)Nmf a 
a 

rfn 



-r 
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where now f*(m, T) is such that /3<&(m, /) = j3mf— E(/) 
is minimum. If m is allowed to assume real values, the 
complexity can be estimated from the knowledge of the 
function 0<f>(m,T) = j3mf*(m,T) - E(/*(m,T)). In- 
deed, it is easy to show that 



d(3$(m,T) 



dm 



(3f*(m,T) 

E(to,T) = E(.T(to,T))=to 2 



dim^P&im^)] (3) 
dm 



= to/3/* (to, T) - /3$(m, T) . 

The function E(/) can be reconstructed from the para- 
metric plot of f*(m, T) and E(to, T). 

Moreover, at fixed to < 1, the glass transition is shifted 
towards lower values of the temperature. Indeed, for 
any value of the temperature T below Tk it exists a 
value to* (T) < 1 such that for to < to* the system is 
in the liquid phase. The free energy for T < Tk and 
to < m*(T) can be computed by analytic continuation 
of the free energy of the high temperature liquid. As 
the free energy is always continuous and it is indepen- 
dent of to < to* (T) in the glass phase (being simply the 
value fmin(T) such that E(/ TO j„) = 0), one can com- 
pute the free energy of the glass below Tk simply as 
F glass (T) = <P(m*(T),T)/m*(T). 

The to copies are assumed to be in the same state. 
This means that each atom of a given replica is close to 
an atom of each of the other to — 1 replicas, i.e., the liq- 
uid is made of molecules of to atoms, each belonging to 
a different replica of the original system. In other words 
the atoms of different replicas stay in the same cage. The 
replica method allow us to define and compute the prop- 
erties of the cages in a purely equilibrium framework, in 
spite of the fact that the cages have been defined origi- 
nally in a dynamic framework. The problem is then to 
compute the free energy of a molecular liquid where each 
molecule is made of m atoms. The to atoms are kept 
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close one to each other by a small inter-replica coupling 
that is switched off at the end of the calculation, while 
each atom interacts with all the other atoms of the same 
replica via the original pair potential. This problem can 
be tackled by mean of the HNC integral equations [28] . 



B. HNC free energy 

The traditional HNC approximation can be naturally 
extended to the case where particles have internal degrees 
of freedom and also to the replica approach where we 
have molecules composed by m atoms. 

We will denote by x — {»!,••• ,x m }, x^ € K d the 
coordinate of a molecule in dimension d. The single- 
molecule density is 



N m 
i—l a—l 

and the pair correlation is 

1,7V m m 

p(x)g(x, y)p(y) = IT ~ ^ IT 5 &jb ~ Uj)) 



(4) 



i,j a=l 



6=1 



(5) 

We define also h(x,y) = g(x,y) — 1. The interaction 
potential between two molecules is v(x, y) — J2 a v (\£-a ~ 
yj). 

"The HNC free energy is given by [13, 28] 



P^[p{x),g{x,y)\ = - J dxdy p{x)p(y)[g(x ,y) log g(x,y) 

- g(x, y) + 1 + 0v(x, y)g(x, y)] 
+ f dxp{x)[\ogp{x) { ^-^[h P T , 



n>3 



(6) 



where 

Tr[hp] n = J dxi ■ ■ ■ dx n h(xi,x 2 )p(x2)h(x2,X3)p(x 3 ) 
■ ■■h(x n -i,x n )p(x n )h(x n ,x 1 )p(x 1 ) . 

(7) 

For Hard Spheres the potential term vanishes, 
J dxdy p(x)p(y)g(x,y)v(x,y) = 0, so the reduced 
free energy (3^ will not depend on the temperature 
in all the following equations. Similarly, all the free 
energy functions that we will consider below do not 
depend on the temperature once multiplied by (3. In 
principle we could stick to /? = 1 and slightly simplify 
the formulae. We have preferred to keep explicitly 
P, in order to conform to the standard notation for 
soft spheres (or for hard spheres with an extra potential). 

Differentiation w.r.t g(x, y) leads to the HNC equation: 
\ogg{x, y) + (3v(x, y) = h(x, y) - c(x, y) , (8) 



having defined c(x, y) from 



h(x, y) = c(x, y)+ dz c(x, z)p(z)h(z, y) . (9) 



The free energy (per particle) of the system is given by 
<t>{m,T) = —— min V\p(x),g(x,y)] , 

Nm p(x),g(x,y) (10) 

$(m,T) = m<f>(m,T) , 

and once the latter is known one can get the free energy 
of the states and the complexity using Eq.s (3). 



C. Single molecule density 

The solution of the previous equations for generic m is 
a very complex problem (it is already rather difficult for 
m = 2). Some kind of ansatz is needed to simplify the 
computation, that may become terribly complicated. 

The single molecule density encodes the information 
about the inter-replica coupling that keeps all the replicas 
in the same state. We assume that this arbitrarily small 
coupling has already been switched off, with the main 
effect of building molecules of m atoms vibrating around 
the center of mass X. & M d of the molecule with a certain 
"cage radius" A. The simplest ansatz for p(x) is then [13] 

p(x) =pj dXj{p(x a - X) , J dup(u) = 1 , (11) 



with 



p{u) = 



e 



(12) 



and p = V 1 j dx p{x) the number density of molecules. 
With this choice it is easy to show that 



Jq J dx P( x ) [ lo § P( x ) ~A = lo § P ~ 1 + 

^(1-to) log(27rA) - ^logm+^(l - to) 



(13) 



D. Pair correlation 

As the information about the inter-replica coupling is 
already encoded in p(x), we make the ansatz for g(x,y): 



g(x,y) = Y[g(\x a -yJ) 



(14) 



where g(r) is rotationally invariant because so is the in- 
teraction potential. We also define G(r) = [g(r)) m . Using 
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the ansatz above, it is easy to rewrite the free energy (6) 
as follows: 

^ = ^Y J drMFoMr^iMr) - [F (r)] m 
+ l + m[F (r)] m - 1 F„(r)} 

dx^)[logp( 2 ;)-l]+-^^Tr[M" , 

n>3 

(15) 



where 



r + u — v\ 



F P (\r\) = J dudvp(u)p(v)g(\r + u- v\)[logg(\ 

F v (\r\) = y dudv_p(u)p(v) g(\r + u — v\)(3v(\r + u — v\) 

(16) 

Note that as g(r) and i>(r) are rotationally invariant, so 
are F p {r) and F v (r). If p(u) is given by Eq. (12), one 
gets 



We will focus first on the g(log g — 1) term in Eq. (6). 
The contribution we want to estimate comes from the 
discontinuity of g(r) in r = D. Thus to compute this 
correction the form of g(r) away from the singularity is 
irrelevant and we will use the simplest possible form of 
9(r). 



A. Expansion of F (r) 

First we will discuss the expansion of F (r) in d = 1. 
The simplest possible form of g(r) is 



g(r)=6(r-D)[l + (y-l)e 



-A»(r-D)l . 



(18) 



the amplitude of the jump of g(r) in r = D is given by y. 
Remember that in our notation r G R and r = \r | € M + . 
As the functions Fq and g are even in r, we can write 



F(\r\) = / du 



where f(r) € {g(r), g(r) log g(r), g(r)(3v(r)}. For Hard 
Spheres F v = 0. 



/oo poo 
dr[Fo(r) m -g(r) m } = 2 dr[F (r) 
-oo JO 

/(|zi + u|) (17) Defining 

erf(i) = — = / rfa; , 



(19) 



(20) 



0(i) = -[l + erf(t)] 



III. SMALL CAGE EXPANSION 

The strategy of [13] was to expand the HNC free energy 
in a power series of the cage radius A, assuming that the 
latter is small close to the glass transition. The expansion 
is carried out easily if the pair potential v(r) and the pair 
correlation g(r) are analytic functions of r. However this 
is not the case for Hard Spheres, as g(r) vanishes for 
r < D and has a discontinuity in r = D, so the formulae 
of [13] for the power series expansion of \& cannot be 
applied to our system. In this section, we will work out 
the expansion in the case where the pair correlation g(r) 
has discontinuities. 

It is crucial to realize, that independently from any ap- 
proximation, in the limit A — > 0, the partition function 
becomes (neglecting a trivial factor) the partition func- 
tion of a single atom at an effective temperature given 
by (3 e ff = 13m. In the case of hard spheres, where there 
is no dependence on the temperature, the change in tem- 
perature is irrelevant. 

In [13] it was shown that the first term of the expan- 
sion is proportional to A if g(r) is diffcrentiable. As we 
will see in the following, in the case of hard spheres, the 
presence of a jump in g(r) produces terms 0(y/A) in the 
expansion. In this paper we will focus on these terms ne- 
glecting all the contributions of higher order in \f~A. This 
means that we can neglect all the contributions coming 
from the regions where g(r) is differentiable and concen- 
trate only on what happens around r = D. 



these functions play the role of "smoothed" sign and 9- 
function respectively; note also that the function 0(t) 
goes to as e~* for t — ► — oo. Then 



f 

1 

2 



e 

du , 9(r + u- D) 



and 



F (r) = 6 



1 + erf 



r-D 



a/47tA 

'r-D 



e - 



= e 



r + D 



r-D 



(21) 



<AA 



+ {y-l)e A ^{e-^ {r - D) Q 



r-D- 2Ap 



4A 



(22) 



r + D 



2Ap\ 



As r > we can neglect the terms proportional to 
^— 2 ^=j in Eq. (22), that give a contribution of order 

exp(— D 2 /A) for A — > 0. Defining the reduced variable 
t = (r- D)/VlA: 



g(t)=6(t)[l + (y-l)e 



F Q (t) = 0(i) + (y- \)e-^ 1 e^ 2 0(t + nVA) 



(23) 
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and Eq. (19) becomes 

/>OC 

/ dr[F (r) m - g(r) m ] = 
Jo 



10 

2VA 



dt[F (t) m -g(t) m ] = 2VAQ(A) 



The last integral, with F a (r) given by Eq. (30) is the 
same as in d = 1, so we obtain 

(24) \ J drF (r) m = X -l dr_G{r) + VAYX d (D)Q m , (32) 



If the function Q(A) has a finite limit Q(0) for A — > we 
will have Q{A) = Q(0) + o(l) and the leading correction 
to the free energy is 0(VAQ(Q)). The limit for A — > of 
Q(A) is formally given by 



/oo 
dt [6(t) m - 6(t) m ] = y m Q r 
-oo 



(25) 



where y m = Y is the jump of G(r) = g(r) m in r = D and 
Q m = di [e(t) m -6(t) m ]. It is easy to show that Q m 
is a finite and smooth function of m for m ^ 0, that 



Q m = (l-m)Q + O[(m-l) 2 } , 

/OO 
d*9(t) log9(t) ~ 0.638 
-oo 



(26) 



and that Q m diverges as Q m ~ yJixjAm for m — > 0. 
Finally we get, recalling that G(r) = [g{r)] m , 

\ J drF (r) m = l -J drG{r) + 2^AYQ m . (27) 

In dimension d > 1 we have, recalling that F (r) and 
G(r) are both rotationally invariant, 

J dr[F (r) m -G(r) m ] = n d J°° drr d ~ x [F (r) m -G(r) m ] , 

(28) 

where fl d is the solid angle in d dimension, Q d = 
2 n d / 2 /T(d/2). The function F (r) can be written as 



Mr) 



du 



{^fA^A) d 



g(\ri + u\) 



(29) 



where i is the unit vector e.g. of the first direction in 
R. . For small \/A the u are small too. The function 
g(\ri + u\) is differentiable along the directions orthogonal 
to i. Expanding in series of u M , fi ^ 1, at fixed u\, 
we see that the integration over these variables gives a 
contribution O(A), so we finally get: 



Mr) 



f 

J — ( 



e 4^ 

d«i -7=g{r + ui) + O(A) , (30) 
V47TJ4 



as in the one dimensional case. The function Fo(r) m — 
G(r) m is large only for r — D^ \f~A so at the lowest order 
we can replace r d_1 with D d ~ x in Eq. (28). We get 

J dr [F (r) m -G(r) m ] = ^ d D d ~ 1 J°° dr [F (r) m -G(r) m ] 

(31) 



where S d (D) is the surface of a d-dimensional sphere of 
radius D, T, d (D) = n d D d ~ l . This result can be formally 
written as 



F (r) m ~ G(r) + 2^AYQ m 8(\r\ - D) 
= G(r)+Q (r) 



(33) 



as the correction comes only from the region close to the 
singularity of g(r), r — D ~ \f~A. 



B. GlogG-term 

Let us now estimate the correction coming from the 
term J drmFo(r) m ~ 1 Fi(r). Using the same argument as 
in the previous subsection, we will restrict to d = 1 . Note 
first that F (r), for \r — D\ ~ \fA, has the form 

F (r)=ye( V -^Jp\ +o(VA) , (34) 



where y is the jump of the function g(r) in r = D. Simi- 
larly, Fi (r) will have the form 



Fi(r) = 



' g{r)\ogg{r) + 0{A) , 



D\ > VA 



,logye(^§) +o(VA) , |r-£>|~^. 

(35) 



The integral 

drlmFoW™- 1 ^^) - m 3 (r) m logg(r)] (36) 



has then two contributions: the first comes from the re- 
gion |r — D\ 3> VA and is of order A as if the function 
g(r) were continuous. The other comes from the region 
\r — D\ ~ \f~A and is of order \f~A as in the previous 
case. To estimate the latter we can use again the re- 
duced variable t and approximate F\(t) <~ y logy @(t), 
F (t) ~y@{t). Then we get 



drimFoir)™- 1 F x (r) - m 3 (r) m logff(r)] - 

y log y 2 VlQ m + o( VI) , 



in d = 1 and finally, in any dimension d, 
i J drmFoirr^FAr) = 



X - J drG(r) log G(r) + VAY log yS d (£>)Q r , 



(38) 
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C. Interaction term 

Substituting Eq. (11) in the last term of the HNC free 
energy one obtains 



Z)C(Z -Y) = H(X -Y)- C(X — Y_), we get 



n>3 



2 ^ n 

n>3 



(43) 



Tr[hp} n =p n J dX 1 ---dX n j 



du\ ■ ■ ■ du n x 



x p(m) ■ ■ ■ piu^hix.! - x_ 2 ,m - u 2 ) 

■ ■ -HX-n - X_ x ,U n - Ml) , 



(39) 



where we used the notations h(X, u) = J|"Li 9(X + u a ) — 
1 and p(u) = nr=i P(Ma) with P(u) g iven b Y E q- ( 12 )- 

The correction 0(y/ r A) to this integral comes from 
the regions where \X, t - X l+l \ = D + 0(/A) for some 
i = 1, • • • , n. In these regions the functions h such that 
their arguments are not close to the singularity can be 
expanded in a power series in u, the correction being 
0(A) [13]. Thus we can write, defining H{r) = G(r) - 1: 

r n Tr[hp) n =JdX 1 ... dX n B(X x - X 2 ) ■ ■ ■ H(X n - X,)- 

n J dX_ 1 ---dX_ n J du\du 2 p(u\)p(u 2 )x 

x [h(X 1 -X 2l ui-u 2 )-H(Xi-X 2 )]x 
xH(X 2 -X 3 )...H(X n -X 1 ) = 

jdX,-.. dX n H{X x - X 2 ) ■ ■ ■ H(X n - X,) 

+ nJdX 1 --- dXM^ - X 2 )x 
xH{X 2 -X^)---H{X n -X 1 ) , 



- NpQ m VAyZ d (D)[H(D) - C(D)} . 
This result is correct in any dimension d. 

IV. FIRST ORDER FREE ENERGY 

Substituting Eq.s (32), (38) and (43) in Eq. (15) one 
obtains the following expression for the HNC free energy 
at first order in \f~A: 

(3F = ^ = 0F o (A) + f3F eq [G(r)} + (3AF[A, G(r)] , 
(3F eq = P -l d d r{G(r) log G(r) - G(r) + 1} 



i r d d k 

2$ J (2^ 
log?- 1 , 



log[l +H(k)]+H(k)--H(kf 



pF = ^(1 - m) log(27r^) + |(1 - m) - ^ logm , 



(3AF = pQ m ^AV d {D) G(D)x 
x [logG(D)-l-H(D)+C(D))] 



where Q m = Qo(l - m) + o((m - l) 2 ), Q 
the Fourier transform has been defined as 



H(k) = p I dre lkr H(r) . 



where in the last step we used Eq. (33): 

du x du 2 p(ui)p(u 2 ) [h(r, u x - u 2 ) - H(r)] = 
F a (r) m - G(r) = Q (r) . 
Collecting all the terms with different n we get 



(40) 



(41) 



(44) 
0.638 and 

(45) 



2 " n 

n>3 



2*-i n 

n>3 



p n TrH n + 



y / dX.dX^XMX, - X 2 )H(X 2 - X 3 )x 
\V(-l) n fP- 3 f dX 4 ---dX n H(X 3 -X 4 )x 



At the first order in \f~A we only need to know the 
function G(r) determined by the optimization of the 
free energy at the zeroth order in y/A, i.e. the usual 
free energy F eq [G(r)\: it satisfies the HNC equation 
log G(r) = H(r) — C(r). Substituting this relation in 
(3AF one simply obtains (3AF = -pQ m VAT, d (D) G(D). 

The derivative w.r.t. A leads to the following expres- 
sion for the cage radius: 

r— 1 — m d , . 

'A* = — — (46) 



Q m pE d (D)G(D) 
which in d = 3 becomes (let us define again Y = G(D)): 
VA* 1 - to 1 



(42) 



D 



Q m 8pY(<p) 



(47) 



x-..H(X n -X 1 ) = 

= 2 £ ^y^/^Tri/" - P - \ dX.dX.d^x 

n>3 n 

x QoiX, - X 2 )H(X 2 - X 3 )C(X 3 - Xj) . 

Substituting the expression of Qo(r) and recalling that 
from the definition of CQT) one has p J dZH(X — 



D 

where ip = - g 9 is the packing fraction. Substituting 
this result in (5AF one has (5AF(A*) = d(m - 1). 

Finally, the expression for the replicated free energy in 
d = 3 is 



/3$(m, <p) = f3F eq (<p) + -(1 - m) log[2^A*(m)] 



+ -(m 



(48) 



■ logm 



nl i I i I i I i I i I i I i fs i I 

tf.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 



8.44 048 052 056 06 064 



FIG. 1: The equilibrium complexity £(<£>) as a function of the 
packing fraction. 



FIG. 2: Phase diagram of the molecular liquid. For m < m* 
(full line) the system is in the liquid phase, for m > m* it is 
in the glass phase. 



Note that for Hard Spheres one has 0F eq (tp) = —S((f), S 
being the total entropy of the liquid. We get then 

dB$ 3 
/?/* (m, <p) = = -| log[27rA*(m)] 



(49) 



^3 . d log 4*(m)^3m-l 

2 dm 2 m ' 

S(m, tp) = m/3f* - p<S> = S(<p) - - log[27rA* (m)] 

3m. .d log A* (to) 3, 



For small enough density the system is in the liquid phase 
and m is equal to 1 at the saddle point. For m = 1 we 
have: 



1 



D 



8Q tpY(tp) 



S mb (tp) = -0t(l,tp) = ^log[2vrA*(l)] 
E(y>) = S(ip)- S vib (tp) 

This allows for a computation of once and 

F(^) are known. Note that \+AtpY(tp) = (3P/p = 
so a model for S(ip) (or V(^)) is enough to determine all 
the quantities of interest. 



V. RESULTS FROM THE HNC FREE ENERGY 

We computed numerically S(tp) and Y(ip) solving the 
classical HNC equation for the Hard Sphere liquid up to 
ip = 0.65. This allows to compute /3$(<p, m) and gives ac- 
cess to all the thermodynamic quantities using Eq.s (49) 
and (50). In this section we discuss the results of this 
computation. We will set the sphere diameter D = 1 in 
the following. 



A. Equilibrium complexity 

The equilibrium complexity is given by Eq. (50). 
It is reported in Fig. 1. We get a complexity S ~ 1 
as found in previous calculations in Lennard- Jones sys- 
tems [9, 13, 26, 27], as well as in the numerical simula- 
tions [9, 10]. The complexity vanishes at ipK = 0.582, 
that is the ideal glass transition density -or Kauzmann 
density- predicted by the HNC equations. 



B. Phase diagram in the [tp,m) plane 

We now compute the thermodynamic properties of the 
glassy phase for <p > tp K . As discussed above, it exists a 
value of m, m*(ip), such that for m < m*(ip) the system 
is in the liquid phase. It is the solution of S(m, tp) = 0, 
where H(m,(p) is given by Eq. (49). In Fig. 2 we report 
m* as a function of tp. Clearly, m* = 1 at tp = px and 
m* < 1 for tp > tpR. m* vanishes linearly at tp c = 0.640. 
As we will see in the following, above this value of tp the 
glassy state does not exist anymore. 



C. Thermodynamic properties of the glass 

The knowledge of the function m*(tp) allows to com- 
pute the entropy of the glass. Indeed, the free energy 
does not depend on m in the whole glassy phase, and 
it is continuous along the line m = m*(tp), so we can 
compute the entropy of the glass simply as 

bgiassKw) = -pFgi ass (tp) = —— (51 

mr{tp) 

This relation is true for m* < 1. Below tp K one has 
m* > 1 and the liquid phase is the stable one. Eq. (51) 
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200 



PP/p 

100 



0.52 0.56 




FIG. 3: Entropy of the liquid (full line) and of the glass 
(dashed line). The two curves intersect at <p k = 0.582 where 
they are tangent and consequently the pressure is continuous 
at the glass transition. The entropy of the glass goes to — oo 
at tp — <p c = 0.640, so the glassy phase does not exist above 
ip c . The dot-dashed line is the entropy of the equilibrium 
states of the liquid, S V ib{<p) = S(p>) — £(<p). 



for to* > 1 gives the entropy of the lowest states in the 
free energy landscape (see below) and can be regarded 
as the analytic continuation of the glass entropy below 
ifK- The reader should notice that the glass phase for 
to* > 1 does not have a simple physical meaning and the 
interesting part of the curves for the glass is in the region 
p > (p K . 

In Fig. 3 we report the entropies of the liquid and the 
glass as functions of the packing fraction. The glass phase 
becomes stable above px — 0.582; note that the entropy 
of the glass is smaller than the entropy of the liquid, i.e. 
its free energy is bigger than the free energy of the liquid. 
The same happens also in Lennard-Jones systems and 
in mean-field spin glass systems. However the physical 
relevant parts of the curves are the liquid one for p < pk 
and the glassy one for ip > pk- 

The reduced pressure, 



pP OS 

— = -f^- 

P oip 



(52) 



is reported in Fig. 4. It is continuous at px and the glass 
transition is a second order transition from the thermo- 
dynamical point of view. Note that the pressure in the 
glass phase is well described by a power law and it has a 
simple pole at p c : 



PPglas 



1 



oc 



(53) 



as one can see from the inset of Fig. 4 where the inverse 
reduced pressure is plotted as a function of p. 

For <p — > p c the pressure of the glass diverges and its 
compressibility x = ^fp vanishes and consequently p c 



FIG. 4: Reduced pressure f3P/ p of the liquid and the glass as 
functions of the packing fraction. The pressure is continuous 
at pk- In the inset, the inverse reduced pressure is plotted; 
in the glass phase it is proportional to <p c — p. 



is the maximum density allowed for a disordered state, 
i.e. it can be identified as the random close packing den- 
sity. The value p c = 0.640 is in very good agreement 
with the values reported in the literature. Note that the 
compressibility jumps downward on increasing tp across 
Pk, i.e. the compressibility of the glass is smaller than 
the compressibility of the liquid. 



D. Cage radius 

The cage radius is given as a function of m in Eq. (47). 
In Fig. 5 we report the cage radius in the liquid phase, 
y/A*(l), see Eq. (50), and the cage radius in the glass 
phase, defined as y/A*(jn*). As Q m ~ y/ir/Am for to ~ 
0, the cage radius vanishes as \J to* for to* ~ 0, i.e. it 
is proportional to y/ip c — f- The vanishing of the cage 
radius for p — > p c means that at p c each sphere is in 



contact with its neighbors, that is consistent with our 
interpretation of tp c as the random close packing density. 



E. Complexity of the metastable states 

From the parametric plot of Pf*(m,p) and S(m, tp) 
given in Eq. (49) by varying to, one can reconstruct the 
function S(/3/) for each value of the packing fraction. 
This function is reported in Fig. 6 for some values of ip 
below and above pk- The function E(/3/) vanishes at 
a certain value {3f m i n , that is given by Eq. (51). The 
saddle-point equation that determines the free energy of 
the equilibrium states is, from Eq. (1), 



dZ(Pf) 
dfif 



1 . 



(54) 
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FIG. 6: Complexity of the metastable states as a function of 
FIG. 5: Cage radius \f~k~ (in units of D) in the liquid and in their free energy /3f for different values of (p. From left to 
the glass phase as function of (p. right, ip = 0.45, 0.5, 0.55, 0.58, 0.6, 0.62, 0.63. The curves are 

truncated arbitrarily at high /3f. The dashed line has slope 1. 



From Fig. 6 we see that this equation has a solution /* > 
fmin for <p < tpx — 0.582. Above ipx Eq. (54) does not 
have a solution so the saddle point is simply /* = f m i n 
and the systems goes in the glass state. In this sense, the 
free energy f m i n of the lowest states below ipx can be 
regarded as the analytic continuation of the free energy 
of the glass, see Fig. 3. The curves £(/?/) in Fig. 6 have 
been truncated arbitrarily at high /3f. We have not done 
consistency checks to investigate where the higher free 
energy states become unstable (i.e. , to compute f m ax)- 



VI. CORRELATION FUNCTIONS 



We will now turn to the study of the pair distribu- 
tion function g(r) in the glass state. In principle a full 
computation would require the evaluation of the correc- 
tions proportional to \[A in the correlation functions of 
a molecule. However we neglect these terms, that we be- 
lieve are small, and we consider again our simple ansatz 
(11), (14) for the correlation function of the molecules, 
in which the information on the shape of the molecule 
is only encoded in the function p(x); these corrections 
should be physically more relevant and interesting. 

As we will see in the following, the correlation function 
of the spheres in the glass is very similar to the one in 
the liquid but develops an additional strong peak (that 
becomes a 5-function at <p c ) around r = D. The integral 
of the latter peak is related to the average coordination 
number of the random close packings. 



A. Expression of g(r) in the glass phase 

We assumed the following form for the pair distribution 
function of the molecular liquid, see Eq.s (11) and (14): 



P-z{x, y) = p(x)g(x, y)p{y) = 

p 2 J dXdY J] p(x a - X)g(\x a - y a \)p(y a - Y) 



(55) 



The pair correlation g(r) of a single replica is obtained 
integrating over the coordinates of all the replicas but 
one: 

9{\xi - 2/i |) = p~ 2 J dx 2 --- dx m dy 2 ■ ■ ■ dy_ m p 2 {x, y) . 

(56) 

Using Eq. (55) we get, with some simple changes of vari- 
able: 



g(r) = g(r) J dudvp(u)p(v)F (\r + u - v|) m_1 , (57) 

where Fo(r) is defined in Eq. (16). The HNC free energy 
is optimized by g(r) = G(r) 1 /™, where G(r) is the HNC 
pair correlation. Thus we get the following expression for 
the pair correlation of a single replica: 



Kr)=G(r) ^ /dM (vS^ Fo(l2: + 2iir " 

u 2 



(58) 



For m = 1, i.e. in the liquid phase, this function is 
trivially equal to G(r). This is not the case in the glass 
phase where m < 1. 



10 



B. Small cage expansion of the correlation function 

We will now expand Eq. (58) for small A. Note first 
that, if r =/= D, the function g(r + u) can be expanded in 
powers of u, and the first correction to g(r) is of order 
A. Then, as before, we will concentrate on what hap- 
pens around r = D. As already discussed in section III, 
around r = D we have, as in Eq. (34), G(r) ~ Y8(r — D) 
and 



F (r) ~y-e 



D 



(59) 



and Eq. (58) becomes 
g(r) = Y9(r - D) / du 



e 

(V4^A) d 



e 



\r + u\-D 



AA 



(60) 

Applying the same argument we used in section III when 
studying the function F (r) in dimension d > 1, we can 
show that the integration over the coordinates u^, ^ ^ 1, 
gives a contribution 0(A). Then we can rewrite, in any 
dimension d: 



g(r) ~ Y6{r - D) 

<-°° dt 



e 4/i 
du^=Q 
VAttA 



u-D 



G(r) 1 



'IT 



AA 

[e(t) 1 "- 1 



m— 1 

■}• 

(61) 



defining the reduced variable f = 



/4A 



The second 



term in the latter expression is a contribution localized 
around r = D. 



C. Number of contacts 

To compute the average number of contacts, let us 
recall that the average number of particles in a shell 
[r.r + dr], if there is a particle in the origin, is given 

by 



dn(r) = Vldr d 1 pg(r)dr . 



(62) 



Thus the number of contacts can be obtained from the 
correlation function g(r). While the full computation of 
the correlation function is rather involved, here we limit 
ourselves to consider the second term in Eq. (61), which 
is proportional to a Gaussian with variance O(VA) that 
becomes a <J(|r| — Z?)-function in the limit A — > 0. 

The value of the number of spheres in contact with the 
sphere in the origin is given by 



(63) 



rD+0(VA) 

z = Q d p rfrr d_1 g(r) 
Jd 



The first term in Eq. (61) gives a contribution 0(y A) 
that can be neglected. If we use r ~ D and G(r) ~ Y at 



the leading order in \/A we obtain, defining the variable 

_ r-D 



z = n d D d - 1 pYx 



r— [°° , f°° dt 

x VAA / de —=e 

Jo J-oo V 77 



-(e-tf 



[e(t) m_1 -i] . 



Recalling that 



1 f 00 

— / dee-^ 2 = Q{t) , 
Jo 



(64) 



(65) 



we get, observing that dt [9(t)— 6{t)\ = 0, and using 
Eq. (46), 



/OO 
dt@{t) [Qit)™- 1 
-OO 

= H d (D)pYV4AQ m = 2d{\ - to) . 



(66) 



This is the expression of the average number of contacts 
at the leading order in v^4, to be computed at to = m* in 
the glass phase. At ip = <p c , where to* = 0, each sphere 
has on average 2d contacts. This is exactly what is found 
in numerical simulations; the condition z > 2d is required 
for the mechanical stability of the packings as can be 
understood by mean of a very simple argument [22] . 

Note that this result is independent on the particular 
expression we chose for S(ip), Y(<p) and G(r), i.e. it 
might hold beyond the choice of HNC equations for the 
molecular liquid provided that the expression (46) for the 
cage radius is correct. 



VII. DISCUSSION 

We will now compare our results with related ones that 
appeared in the literature. The main obstacle for a quan- 
titative comparison is that the HNC equations are known 
to yield a not very good description of the Hard Sphere 
liquid at high density [28]; typically one would obtain 
the right curves if one shifts the value of ip of a quantity 
of order 0.03. Therefore, we should limit ourselves to 
a qualitative comparison of the results coming from the 
HNC equations with the results of numerical simulations. 
However, note that, although the expressions (47), (48) 
for the replicated free energy have been derived start- 
ing from the expression (6) for the HNC free energy, the 
final result depends only on the equilibrium entropy of 
the liquid S(ip). It is interesting then, for the purpose 
of comparing our results with experiments and numer- 
ical simulations, to consider a more accurate model for 
S(ip) in the liquid phase. We repeated the calculations 
of section V substituting the Carnahan-Starling (CS) en- 
tropy [28] 



Scs(<p) = - log 



6</A 4(p -Zip 2 



ire J (1 — ip) 2 



Y CS (<P) = 



(1 - p)3 ' 



(67) 
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FIG. 7: Equilibrium complexity T,(ip) as a function of the 
packing fraction. The full line is from the HNC equation 
of state (see Fig. 1), the dashed line is from the Carnahan- 
Starling equation of state. The black dots are numerical data 
of Angelani et al. [15] (the data shown here correspond to 
ag 1 ' in Ref. [15]), the dot-dashed line is extrapolated from 
the numerical data reported by Speedy [3]. 

instead of the HNC entropy in Eq.s (48), (47). All the re- 
sults of section V are qualitatively reproduced using the 
CS entropy, but the latter gives results in better agree- 
ment with the numerical data. However, this procedure 
is not completely consistent from a theoretical point of 
view: one should always keep in mind that our aim here 
is not to present a quantitative theory, but only to show 
that the replica approach yields a reasonable qualitative 
scenario for the glass transition in Hard Sphere systems. 



A. Complexity of the liquid and Kauzmann density 

In Fig. 7 we report the equilibrium complexity E(<p) 
obtained substituting the HNC and the CS expression 
for S(<p) and Y(ip) in Eq. (50). The results are compared 
with recent numerical results of Angelani et al. [15] ob- 
tained on a 50 : 50 binary mixture of spheres (to avoid 
crystallization) with diameter ratio equal to 1.2: the vi- 
brational entropy was estimated using the procedure de- 
scribed in [9, 29] and the complexity was computed as 
S(ip) — S v ib(if). A quantitative comparison is difficult 
here because in the case of a mixture there can be correc- 
tions related to the mixing entropy, S m i X ~ log 2. Never- 
theless the data are in good agreement with our results. 
A detailed comparison would require the extension of our 
computation to binary mixtures following [9]. 

Another numerical estimate of E(y>) was previously re- 
ported by Speedy [3] , who rationalized his numerical data 
assuming a Gaussian distribution of states and a partic- 
ular form for the vibrational entropy inside a state. The 
free parameters were then fitted from the liquid equation 
of state. The curve obtained by Speedy also agrees with 



0.12 




FIG. 8: Inverse reduced pressure -^p of the Hard Sphere liquid 
as a function of tp. The black dots are from the simulation of 
Rintoul and Torquato [1]. The full line is obtained from the 
CS equation of state while the dot-dashed line is from the 
HNC equation of state. The dashed parts of the two curves 
correspond to the (ideal) glass phase. Note that all the curves 
are quasi-linear functions of ip in the glass phase. 



our results. 

Both the HNC and the CS estimates of the Kauzmann 
density (ifK = 0.582 and ipx = 0.617 respectively) fall, 
as it should be, between the Mode-Coupling dynamical 
transition that is <pmct ~ 0.56 [5, 6], and the Random 
Close Packing density that is estimated in the range <p = 
0.64-^0.67, see e.g. [16]. 

A computation of E(y>) based on very similar ideas 
was presented in [14], where a very similar estimate of 
ifiK ~ 0.62 was obtained. However in [14] the complexity 
was found to be E ~ 0.01, i.e. two orders of magnitude 
smaller than the one obtained from the numerical sim- 
ulations. This negative result is probably due to some 
technical problem in the assumptions of [14]. 



B. Equation of state of the glass 

In Fig. 8 we report as black dots the numerical data for 
the pressure of the Hard Sphere liquid at high ip obtained 
by Rintoul and Torquato [1]. The data were obtained 
extrapolating at long times the relaxation of the pressure 
as a function of time after an increase of density starting 
from an equilibrated configuration at lower density. We 
also report the curves of the pressure as a function of the 
density obtained from the HNC and CS equations, both 
in the liquid and in the glass state. 

The agreement of the HNC curve with the data is not 
very good even in the liquid phase, due to the modest 
accuracy of the HNC equation of state. However, the 
qualitative behavior of our curve is in good agreement 
with the numerical data, and in particular the quasi- 
linear behavior of the inverse reduced pressure in the 
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glass phase found in [1, 3], -Sp oc tp c — Lp, is reproduced by 
the HNC curve. The HNC pressure of the glass diverges 
at ip c — 0.640 as discussed in section V; the latter is the 
HNC estimate of the random close packing density. 

The CS curve describes well the pressure in the liquid 
phase [28] . Comparing the curve with the data of Rintoul 
and Torquato, we see that the glass transition happens in 
the numerical simulation at a density ip g ~ 0.56 smaller 
than the one predicted by the CS curve, ipK = 0.617 [30], 
and very close to the Mode-Coupling transition density, 
fMCT ~ 0.56. This is not surprising, since the relaxation 
time grows fast on approaching the ideal glass transition; 
at some point it becomes larger than the experimental 
time scale and the liquid falls out of equilibrium becoming 
a real glass. It is likely that the data of Ref. [1] describe 
the pressure of a real nonequilibrium glass, while our com- 
putation gives the pressure of the ideal equilibrium glass, 
that cannot be reached experimentally in finite time. 

C. Random close packing 

Both the HNC and CS equations predict the existence 
of a random close packing density ip c where the pressure 
and the value of the radial distribution function g(r) in 
r = D diverge. The HNC estimate is ip c = 0.640, in the 
range of the values (ip c — 0.64 0.67) reported in the 
literature. The CS estimate is <p c = 0.683 and it is also 
a value consistent with numerical simulations. 

The reader should notice that the theoretical value for 
(p c is related to the ideal random close packing; how- 
ever the states corresponding to this value of <p c can be 
reached by local algorithms, like most of the algorithms 
that were used in the literature, in a time that should 
diverge exponentially with the volume. Some caution 
should be taken in using the data obtained by numerical 
simulations. The question of which is the value of the 
density that can be obtained in large, but finite amount 
of time per particle is very interesting and more relevant 
from a practical point of view: however we plan to study 
it at a later time. 

Note that the computation of the mean coordination 
number z of section VI, that gives z = 6 at <p — ip c in 
d = 3, is independent of the particular form we choose 
for S(ip), and thus is valid for both the HNC and CS 
equations of state. The value z = 6 has been reported in 
many studies [19-23]. 

VIII. CONCLUSIONS 

We successfully applied the replica method of [13, 25] 
to the study of the ideal glass transition of Hard Spheres, 



and in general of potentials such that the pair distribu- 
tion function g(r) shows discontinuities, starting from the 
replicated HNC free energy and expanding it at first or- 
der in the cage radius \J~A. 

This result allowed us to compute from first princi- 
ples the configurational entropy of the liquid as well as 
the thermodynamic properties of the glass up to the ran- 
dom close packing density. Our computation is based 
on the HNC equation of state, that is known to yield a 
poor quantitative description of the liquid state at high 
density. Nevertheless, we found that the qualitative sce- 
nario for the ideal glass transition that emerges from the 
replicated HNC free energy is very reasonable. In par- 
ticular, we found a complexity S ~ 1, a Kauzmann den- 
sity tpx — 0.582, and a random close packing density 
(p c = 0.64. All these results compare well with numerical 
simulations. 

Using, on a phenomenological ground, the Carnahan- 
Starling equation of state instead of the HNC equation of 
state as input for our calculations, we could also compare 
our results with the high-density pressure data of Rintoul 
and Torquato showing that they are indeed compatible 
with the observation of a real glass transition. 

Moreover, we found that the mean coordination num- 
ber in the amorphous packed states is z = 2d irrespective 
of the equation of state we use for the liquid, in very good 
agreement with the result of numerical simulations and 
with theoretical arguments [19, 20, 22, 23]. 

It is worth to note that our results do not prove the 
existence of a glass transition for the Hard Sphere liquid, 
as they derive from a particular approximation for the 
molecular liquid free energy (the HNC approximation), 
and, in general, other approximation such as the Percus- 
Yevick are possible [28]. 
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